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Abstract 



Abrikosov's solution of the linearized Ginzburg-Landau theory describing a periodic 
lattice of vortex lines in type-II superconductors at large inductions, is generalized to 
non-periodic vortex arrangements, e.g., to lattices with a vacancy surrounded by relaxing 
vortices and to periodically distorted lattices that are needed in the nonlocal theory of 
elasticity of the vortex lattice. Generalizations to lower magnetic inductions and to three- 
dimensional arrangements of curved vortex lines are also given. Finally, it is shown how 
the periodic vortex lattice can be computed for bulk superconductors and for thick and 
thin films in a perpendicular field for all inductions B and Ginzburg-Landau parameters k. 



1 Introduction 

From Ginzburg-Landau (GL) theory [TJ Landau's thesis student Alexei Abrikosov predicted 
that superconductors with a GL parameter k > l/y/2 may contain a lattice of vortices of 
supercurrent, or flux lines (fluxons) with quantized magnetic flux $o = h/2e = 2 • 10~ 15 Tm 2 . 
Abrikosov had linearized the GL equations with respect to a small order parameter \tp\ 2 and 
discovered a solution ip(x,y) possessing a regular lattice of zero lines. This lattice solution 
appears when the applied magnetic field B a (along z) is decreased below the upper critical field 
B c2 = $o/(27r£ 2 ), where £ = X/k is the GL coherence length. At B a = B c2 one has the average 
induction B = B c2 . With decreasing B a , the induction decreases and reaches B = at the lower 
critical field B a = B c \ = $o(ln/t + a)/(47rA 2 ) with a(n) rs 0.5 for k ^> 1 (see below). At the same 
time the vortex lattice spacing a m (Qq/B) 1 / 2 increases and diverges when B — > 0. Abrikosov 
tells that he had obtained this vortex solution in 1953 but Landau didn't like it, stating that 
there are no line-like singularities in electrodynamics. Only when Feynman [2] had published 
his paper on vortices in superfluid Helium, did Landau agree and Abrikosov could publish his 
solution in 1957 [3]. For his prediction of the vortex lattice Abrikosov 50 years later in 2003 
received the Nobel Prize in Physics together with Vitalii Ginzburg and Anthony Leggett. 

After first evidence of the triangular vortex lattice in superconducting Niobium by small-angle 
neutron scattering in Saclay [I], Trauble and Essmann in Stuttgart succeeded [5j[6] to observe the 
vortex lattice directly by decorating the surface of a superconductor with iron microcrystallites 
("magnetic smoke"). At that time I joined this research group headed by A. Seeger and wrote 
my thesis on the theory of defects in the vortex lattice [Zj. Parts 1 and 2 deal with low inductions 
B <C B C 2 when London theory may be used and the vortices interact with each other pairwise, 
similar to 2D atomic lattices. Parts 3 and 4 consider high inductions B « B 2 , where the shape 
of the GL solutions ip(x,y) and B(x,y) may be obtained from linearized GL theory, while the 
nonlinear GL terms determine the amplitudes of this ip and B. My thesis extended Abrikosov's 
theory of periodic vortex lattices to non-periodic vortex arrangements, see below. Such distorted- 
lattice solutions are required to calculate the elastic energy of the vortex lattice and the energy 
of lattice defects like vacancies and dislocations. They are also helpful to visualize where the 
solutions of the linearized GL theory apply and how they have to be modified at lower inductions. 
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2 Abrikosov's ideal vortex lattice near B C 2 



In the usual reduced units (length A, induction \/2B c , energy density B 2 /fx , where B c = 
B c2 /y/2n is the thermodynamic critical field) the spatially averaged free energy density F of 
the GL theory referred to the Meissner state (ip = 1, B = 0) within the superconductor reads 
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Here ^(r) = / exp(iyj) is the complex GL function, B(r) = V x A the magnetic induction, 
A(r) the vector potential, and (...) = (1/V) J v d 3 r . . . means spatial averaging over the super- 
conductor with volume V. Introducing the super velocity Q(r) = A — Vcp/n and the magnitude 
/(r) = \ifj\ one may write F as a functional of the real and gauge-invariant functions / or f 2 = u 
and Q, 

F = (^P^ + + fQ 2 + (VxQ)j . (2) 

In the presence of vortices Q(r) has to be chosen such that VxQ has the appropriate singularities 
along the vortex cores, where / vanishes. By minimizing this F with respect to ip, A or /, Q, 
one obtains the GL equations together with the appropriate boundary conditions. For the 
superconducting film considered in Sec. 5, one has to add the energy of the magnetic stray field 
outside the film, which makes the perpendicular component B z of B continuous at the film 
surface, see below. 

The two GL equations are obtained by minimization of F (1) with respect to ip and A, 
5F/5ip = and 5F/5A = 0, yielding 

(VA-/tA) 2 ^ = /t(l-|^| 2 )V, (3) 
V x [V x A] = H 2 Q. (4) 

With B a and B chosen along the z axis and in the gauge A x = —By + A x (x, y), A y = A y (x,y) 
(A x , A y are terms of higher order) the linearized first GL equation, obtained by omitting the 
term \ijj\ 2 in (3), has the general solution 

ij>(x, V) = exp(-nBy 2 /2) g(x, y) , (5) 

— + i— = (6) 

dx dy 

This means g(x,y) = g(z), z = x + iy, can be any analytical function. For a periodic solution 
satisfying |^| 2 (r + R mn ) = |^| 2 (r), r = (x,y), with real and reciprocal lattice vectors 

R-mn = {mxi + nx 2 ; ny 2 ) , (7) 
K mn = (2ir/x 1 y 2 ){my 2 ; -mx 2 + nx x ) , (8) 

(m,n = 0, ±1,±2, triangular lattice: X\ = a, x 2 = Xi/2, y 2 = Xiv3/2; square lattice: 
Xi = y 2 = a, x 2 = 0) and with a zero at r = 0, one obtains for g(z) the function i?i defined as [8] 

oo 

$i{z,t) = 2^(-) n exp[i7rr(n + |) 2 ]sin(2n + l)z . (9) 

n=0 

Thus, the periodic Abrikosov solution with zeros at the r = R mn may be written as 

Mx,y) = exp (- ^(—(x + iy), X2 + m ) . (10) 

v xiy 2 ^ v xi x x ' 
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This solution has the mean induction B = $0/(^12/2), normalized order parameter (IV^I 2 ) = 1, 
and the Fourier series \4>a\ 2 = u A (x, y), 

Ua (t) = Yl {-) mn+m+n exp ( m^l^ e lKr . (11) 

mn 

From the zero u>a{0, 0) = follows that the sum over all Fourier coefficients in (11) is zero for 
all lattice symmetries. (Abrikosov [3] chose a different position for uja = and thus obtained 
the function $3 [8]). 



3 Distorted vortex lattice near B C 2 

The GL solution for a distorted vortex lattice near B c2 is obtained as follows. Assume that each 
of the straight and parallel vortex lines is displaced from its ideal lattice positions R mn = R„ = 
(X V ,Y V ) by displacements s„ = (s ux ,s uy ), r u = (x v ,y v ) = H u + s u , such that the displacement 
field itself is periodic with a super lattice N times larger than the vortex lattice, but with same 
symmetry, s(r + NH mn ) = s(r). Where needed we use a continuous displacement field s(r) 
defined such that it has the same Fourier transform as the discrete s u = s(R iy ). A distorted 
triangular vortex lattice with spacing x\ = a then has the solution, Eq. (5), 

«(«, V) = c, exp (- %) r) n M ;! l N f N iZ - Z "-^Y ] < 12 > 

v V3a 2/ a V i?i[(7r/iVa)(2 - z v ),t\ 

with z = x + iy, z u = X u + iY u , r = (1 + i^/3)/2, s u = s xv + is yV) the product is over one super 
cell, and c\ ~ 1 is a normalization constant. When all s v = 0, the product in (12) is unity, 
il^ = 1, thus the first two factors in (12) are the ideal lattice solution with c\ = 1, cf. Eq. (10). 
Each factor of the product shifts one zero from r = R y to r = R^ + s^. The absolute value of 
\ip\ 2 = io of this GL function may also be expressed in terms of the Fourier series ^(r) (11), 

" (r| = C '" j(r| n ^[(r-R„)/iV] ' (13) 

In the limit of infinite super cell, N — > 00, one may use $i(z/N,t) oc z/N for \z\/N 1, thus 
one may replace the function $i by its argument since all the constant factors cancel or combine 
to a normalization factor that follows from numerics. One then obtains simply 

^)=c\ 4)n '7 R ;"i" 12 ( i4 ) 

r - R„ r 



4 Vortex lattice vacancy near B, 



c2 



Removing the central vortex at R„ = adds a factor 1/r 2 to the linearized solution w(r). 
Obviously, if the other vortices are not allowed to relax, this solution at large distances vanishes 
as 1/r 2 ; it cannot be normalized and its energy is infinite. However, if the relaxation of the other 
vortices is chosen appropriately it will minimize the defect energy and make it finite. This can 
be seen from the solution 



oj r 



uj A (r) h(z) 
h(0) 
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(15) 
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Figure 1: Order parameter u(x, 0) for the ideal vortex lattice (dashed line) and for vortex lattice 
with vacancy, Eq. (15), with simple relaxation field (19) (dotted line) and with better relaxation 
field that minimizes the defect energy (solid line, see text). 



/5 = S4 = ^f> 1 - (16) 



The constant factor 1 1 / h(0) | 2 was added to force convergence of the infinite product. The solution 
for a super lattice of vacancies positioned at the NH U , is given by expression (13) divided by 
u(t/N) that removes the zeros at positions NH U .. 

The energy of both the ideally periodic and the distorted vortex lattices is calculated via the 
Abrikosov parameter 

This (3 enters the free energy of the linearized GL theory (referred to the normal state), that has 
to be minimized when B is held constant, 

F _B 2 (B c2 - BY 

2/xo 2/i [l + (2«: 2 -l)/?] 1 ' 

and the free enthalpy that has to be minimized when B a is held constant, 

r = p BB a = (B c2 - B a f 

/i 2// (2k 2 - 1)(3 ■ 1 ] 

The elastic energy of the distorted vortex lattice is the product of the derivatives OF/ 8(3 or 
dG/d(3 times the change of f3{ip} times the volume, with the limit of infinite volume taken. This 
means that all elastic energies and energies of structural defects near B c2 vanish as (B c2 — B) 2 oc 
(B c2 — B a ) 2 . This is true also for the shear modulus Cqq of the vortex lattice, which can be 
obtained using Abrikosov's periodic lattice solution [7J [9] . 

In the case of the vortex vacancy, the resulting defect energy is finite only if the vortices relax 
(shift towards the removed vortex) such that at large distances (and after numerical minimization 
practically at all distances) the vortex displacements are 

^ - -ass < 19 > 
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Figure 2: Contour lines of the order parameter u(x, y) (15) of a vortex lattice with one vacancy 
at x = y = and complete relaxation, see solid line in Fig. 1. The vortex displacements are 
indicated by short bold lines between two dots. 

with n = -B/$o = 1/(^12/2) the vortex density. If the radial displacements were chosen smaller 
(larger) than in (19), the order parameter (15) would vanish (diverge) at large distances r. But 
with the correct displacements that minimize f3 and thus the defect energy, the amplitude of the 
oscillating order parameter stays almost constant, even near the vacancy. This can be seen in 
Fig. 1, where the profiles of u(x,y) along y = are plotted for ideal triangular vortex lattice 
and for the lattice with a central vacancy with the simple relaxation (19) and with an improved 
relaxation field = — H^y/Ea 2 /(AttR 2 ,) + 0.068a 4 / Rf]. Figure 2 shows the contour lines of the 
fully relaxed order parameter u(x,y), which has a maximum at the origin (vacancy position) 
and minima (zeros) at the vortex positions. 

The displacement (19) means that the area of the "relaxing ring", 2ttR i/ s u = 1/n = x\yi = 
<&o/B, exactly equals the area of one lattice cell. In other words, the removal of one vortex 
at Rj, = is compensated by the relaxation of the surrounding vortices such that the average 
vortex density in any contour (containing or not containing the vacancy) stays constant and 
equals the density n that was there before the vacancy was introduced. Note that the field (19) 
satisfies V • s(r) = 0, and thus describes a pure shear deformation. More precisely, one has 
V ■ s(r) = (1/71)02(1") (02 is the 2D delta function), i.e., the displacement field (19) "remembers" 
that one vortex cell area was removed. Further interesting properties of structural defects in 
vortex lattices and other two- or three-dimensional soft lattices are discussed in [10]. 

5 Distorted vortex lattice away from B c2 

As shown with the vacancy example, the distorted-lattice solution (14) of the linearized GL 
equations yields finite energies of lattice defects only if most of the vortex displacements are 
allowed to relax appropriately. But without this relaxation, the elastic energy is infinite. For 
example, if only the one vortex at the origin is displaced, s u = so^ox (^0 = 1 if ^ = 0, else 
5 v q = 0) one has from (14) 

u{r) = cu A (r)\r - s\ 2 /r 2 = u A ■ (1 - 2xs /r 2 ) + 0{s 2 ) , (20) 
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i.e., the periodic order parameter is modulated by a slowly decreasing function. The Abrikosov 
(3 of this defect times the volume, limy^ 00 ( ( 5 — /3 )V, diverges and so the defect energy diverges. 

This unphysical divergence of defect energies of the vortex lattice is removed when the in- 
fluence of the nonlinear GL terms on the solutions u>(x,y) and B(x,y) is accounted for. This 
calculation was performed in a series of 4 papers [11J: Parts 1 and 2 deal with the linear elastic 
energy of the vortex lattice at low and high inductions B. Parts 3 and 4 derive the GL solu- 
tions for the distorted vortex lattice when the vortex lines are straight and parallel or arbitrarily 
curved. The essential result is that the long-ranging modulation factors like (1 — 2s/r) in the lin- 



earized solution (14) become exponentially damped over a new length £' = 1/k^ = £/y2(l — b) 
with b = B/B C 2. As b — > 1, this screening length becomes infinite and the linearized solution 
(14) is recovered. At b < 1, the distorted-lattice solution (14) should be replaced by 



2 



u(r) =u A (r)[l + J2^VK (\r-R u \k i ,)\ + 0(s 2 ), (21) 

where Kq(x) is a modified Bessel function with the limits Kq(x) ~ — lnx (x <C 1), Kq(x) ~ 
(n/2x) 1 ^ 2 e~ x (x ^> 1). This generalized expression up to terms linear in the vortex shifts 
reproduces the linearized solution (14), (20) when kf — ► 0, but it does not possess the correct 
zeros at r v = + s„. This may be corrected by replacing in (21) the periodic order parameter 
uja(t) by the "phase modulated" ua[t — s(r)] and cutting the infinity of K off. The resulting 
solution is still exact up to linear terms in s u since Vu>A(r) vanishes at the and thus the 
expansion of (jJ A \? — s(r)] contains no linear term. 

The screening length £' = 1/k^ may be derived by considering only one Fourier component 
of the displacement field, 

s v = Re{s exp(ikR„)} (22) 
with k = (k x , k y , 0) and Re = real part. One may then write the linearized solution as 

w(r)= Wy4 (r)[l + ir,(r)] 2 + 0( S 2 ), (23) 

-iW = 2 E = '^w exp|!(k + K ,r| } ■ (24) 

In r](r) (24) the terms with reciprocal lattice vectors K / shift the zeros of w(r) ("phase 
modulation"), while the term K = yields an "amplitude modulation" of lu(t). This term 
diverges as 1/k 2 , i.e., it yields a diverging amplitude modulation when the wavelength of the 
displacement field is large. 

From physical reasons it is clear that this term oc 1/k 2 has to be cut off, e.g., replaced by 
1/ (k 2 + k 2 ^). Accounting for all the GL terms nonlinear in u oc 1 — b (b = B/B C 2, terms like u 2 , 
B 2 , Q 2 ) indeed yields such a cut off, with k 2 ^ = 2(1 — b)/^ 2 . The resulting solution for periodic 
s(r) may be written as 



2b Vs(r 

1 + 



£ 2 k 2 + k 2 



+ 0(s 2 ). (25) 



co(r) = uj A [r - s(r)] 

In a similar way, the solution for the induction B(x,y) of the linearized GL theory 

(cu)_^w(r) 
2k 



B(t) = B + B c2 ^^> (26) 



is modified by the nonlinear terms to give for periodic s(r) 



B(r) = B [r-s(r)]-^^ + O(s 2 ) (27) 
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with k\ = 1/X' 2 = (lo}/\ 2 ~ (1 — b)/X 2 and B (x,y) the ideal periodic solution for s = 0. In 
deriving (27) all terms containing k^ have cancelled. From the solutions (25) and (27) for periodic 
s(r), the generalization to arbitrary displacement fields is obtained by Fourier transform. 



6 Curved vortices 

The above method can be extended to 3D displacement fields s v (z) = [s vx (z), s uy (z), 0] describing 
distorted lattices of curved vortices, 

s "( z ) = I F~T" § ( k ) ex P(^ kR ^) , 

jbz 87r d n 

s(k) = ldzs u {z) exp(-ikR„) , (28) 

where now r = (x, y, z), R„ = (x u , y v , z), n = B/§ , and the k integration extends over the first 
Brillouin zone of the ideal vortex lattice [since s(k + K) = s(k)] and over — oo < k z < oo. The 
coordinate z plays here the role of a line parameter. The order parameter which solves the GL 
equations near B c2 and has zeros at the vortex positions r v (z) = H u + s v (z) is 



uir) = u a [t) 



1 + 



E fdz's v (z') 



exp(- 


|r-R'JM 


2\v-K 





+ 0{s 2 ). (29) 



The 3D solution for the induction B(r) = V x A(r) for periodic s(r) after averaging over a 
vortex cell may be written as 

B(,) = ifi + g ' V ^ + + ff/ fe +0 (^), (30) 

A(r) = ij3zxr + B 1 !^A + 0( S 2 ). (31) 

These expressions coincide with the first-order expansion terms (in s) of a linear superposition 
of spherical "source fields" centered at each vortex element: 

B(r) = « o4 !E/*^^, (32) 

where p m [(r — r u ) 2 + a 2 /4] 1//2 has an inner cut-off a/2, half the vortex spacing in our 
derivation from the vortex lattice. The expression (32) is also the solution of London theory for 
arbitrarily arranged curved or straight vortices if one puts kh — 1/X (i.e. 6^0) and the vortex 
core radius r c m £ for the inner cutoff. The line element of the path integral in (32) may be 
parameterized with z as line parameter and integration variable, 

dr v =*^dz=(*+*p.)dz. (33) 
dz v dz ' 

7 Nonlocal elasticity of the vortex lattice 

The distorted-lattice solution up to terms linear in the displacements s can be used to calculate 
the linear elastic energy of the vortex lattice, F elast = -F{s^} — F{s u = 0}, referred to the perfect 
lattice (the equilibrium state). The most general expression quadratic in the 2D displacements 
s v (z), or in their Fourier transforms s(k) = (s x ,s y ,0), (28), is 

1 r d 3 k 

i^ciast = ^ y BZ s a (k)$ Q/3 (k)5 /3 (-k) , (34) 
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where the sum over the indices a, (5 — (x,y) is taken. The 2x2 matrix $ a p is the elastic 
matrix. This expression applies for both an elastic continuum and for a lattice. For a lattice 
§ a/ 3 is periodic, $ Q/ 3(k + K) = $ ai g(k), and thus the integral should be restricted to the first 
Brillouin zone (BZ). The BZ for the triangular lattice is a hexagon, and for the square lattice a 
square. Where required, the BZ may be approximated by a circle with radius k B = (26) 1 / 2 /^, 
b = B/B c2) and area nk^ = 4ir 2 n, n = B/Q . 

For a uniaxial elastic continuum the elastic matrix Q a p{k X) k y , k z ) has the form 

n$ a/3 (k) = (c n - c && )k a kp + 5 afS [(k 2 x + k 2 y )c && + k 2 z c 44 ] . (35) 

In it the coefficients are the elastic moduli: Cn — c 66 the isotropic compression modulus, Cn 
the uniaxial compression modulus, Cqq the shear modulus, and C44 the tilt modulus. The elastic 
moduli of the vortex lattice are obtained by deriving the elastic energy, e.g., from GL theory 
and comparing it at k 2 + k 2 <^. k 2 B with the continuum limit (35). This yields 

_ B 2 dB a 1 

Cll{k) ~ ^W(iT«MW + c " (36) 

BB c2 2 (2re 2 -l)2re 2 

C66 = 8^ (1 ~ 5) ^TTT7W ( } (37) 

<*(*) = -T-^7p+^^- (38) 
/i l + k 2 /k 2 h no 

These expressions are exact at large reduced induction b = B/B c2 — > 1 and for all re, but they 
are written such that they reduce to the correct values also in the limit of small induction 
B <C B c2 . In cqq, (3a = 1.160 is the Abrikosov parameter of the triangular lattice (the square 
lattice is unstable and thus has negative c 66 ); the third factor reduces to 1 for 2re 2 1 and to 
(2re 2 — 1)0 A — > for re — > l/y/2, which means the shear stiffness of the vortex lattice is zero in 
superconductors with re = 0.71; the factor 1 — 0.36 interpolates between the the correct limits 
at b — > 1 and b — > 0. In particular, for b <C 1 and 2re 2 3> 1, (37) reproduces the London result 
c 66 = BB c2 /(8k 2 /jLo). 

An interesting result is the dependence of cn (36) and C44 (38) on fc = |k|, which means 
the elasticity of the vortex lattice is non-local. In the limit of uniform stress, k — > 0, these 
expressions reproduce the known values of the compression and tilt moduli obtained by thermo- 
dynamics, Cn — c 66 = [B 2 1 fj, )dBa/dB, c 44 = BB a /fj, . However, when the wavelength of the 
periodic compression or tilt decreases, i.e., the wave vector k increases, these moduli decrease. 
This means, the vortex lattice is softer for short- wavelengths compression and tilt than it is for 
long wavelengths. The two characteristic lengths or wave vectors were already introduced above, 



k h = 1/A' w v / l 3 6/A and k^ = 1/? = ^2(1 - &)/£ . 

This dispersion or elastic non-locality means, e.g., that a point force exerted by a small 
pinning center on the vortex lattice, deforms the vortex on which it acts not like plugging a string 
but more, causing a sharp cusp since a local deformation costs little energy. If the interaction of 
the vortices with the pinning center is via the order parameter \ip\ 2 or via the gradient term in 
the GL functional, then this interaction itself is nonlocal, smeared over the length £' = \ jk^. In 
the expressions for the elastic force and the elastic energy there is thus a factor 1 + k 2 jk\ in the 
numerator that compensates the same factor in the denominator originating from cn(k), (36). 
Therefore, the factor 1/(1 + A; 2 /A; 2 ) in cn has no physical meaning in pinning problems since 
near B c2 where £' can be larger than the vortex spacing a, it is not possible to exert a pinning 
force on one single zero of the order parameter but only on an area with radius £' containing 
several such zeros. The nonlocality factor 1/(1 + k 2 jk\) in C44, however, is important in pinning 
theories since it strongly enhances the elastic deformations caused by small pins acting on the 
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vortex cores. In (not very realistic) models where the pinning force acts only the magnetic field 
of the vortex but not on the vortex cores, this enhancement of the elastic displacement may 
vanish, cancelled by the non-locality of this model force. 

The correct, non-local elasticity thus effectively softens the vortex lattice and leads to large, 
pinning-caused distortions and disorder of the vortex lattice. Furthermore, the thermal fluctu- 
ations of the vortex lattice are strongly enhanced by this non-local elasticity. In both cases the 
lattice softening is caused mainly by the dispersion of 044(h), while the dispersion and reduction 
of C\\(h) is not so important since the shear modulus Cqq is typically much smaller than Cu(h) 
and the shear modes of the elastic deformation thus dominate over the compressional modes. 
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Figure 3: The magnetization curves of the triangular vortex lattice (solid lines, numerical result), 
coinciding within line thickness with those of the square lattice. Shown are h = B a /B c2 versus 
b = B/B c2 (upper left triangle) and — m = h — b versus h (lower right triangle). The dots show 
the fit, Eq. (59), good for k < 20. 



8 Vortex arrangements at low inductions 

At low inductions b < 0.2 and not too small k > 2, the GL theory for arbitrary 3D arrangements 
of vortices reduces to the London theory, which may be expressed by the energy functional 

F{B} = Jd 3 r[B 2 + A 2 (V x B) 2 ] . (39) 

Here A is the London depth equal to the GL magnetic penetration depth. Minimizing -F{B} 
with respect to the induction B(r) using VB = 0, and adding appropriate singularities along 
the positions v v (z) of the vortex cores, one obtains the modified London equation [12J, 

(_A 2 V 2 + l)B(r) = $0 E / d ^ 5 3(r - r„) (40) 
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Figure 4: The magnetic field B(r) and order parameter IV^ 7 ")! 2 of an isolated vortex line calcu- 
lated from Ginzburg-Landau theory for GL parameters k = 2, 5, and 20. For such large k the 
field in the vortex center is twice the applied equilibrium field, -B(O) « 2B c i = 2B a . 



with 5 3 the 3D delta function. From this one obtains the energy of an arbitrary arrangement of 
straight or curved vortices, 

fW«)} = 3-t|— EE faJ^ gHgA) . (4i) 



In this double sum the terms describe the pairwise interaction of the vortex line elements 

G?r M , dr„ over the distance = |r M — r„|. The term = v is the self-energy of the /ith vortex 
line, which depends on the shape of this vortex. In it an inner cut off is needed, obtained, e.g., by 
putting r^Jz, z') = |r M (z) — r^(z')\ 2 + r\ with r c « £ the vortex core radius, to avoid divergence 
when in the integral the parameters equal, z = z'. 

From the GL nonlocal elastic energy (34)- (38) one may construct an effective interaction 
potential between vortex line elements such that the full nonlocal linear elastic energy is repro- 
duced at small displacements [13]. At the same time, this interaction at low b <^ 1 reproduces 
the London interaction for arbitrary vortex arrangements, and an approximate GL interaction 
valid at all b and k, 



2 



exp(-?>/A') r f exp(-7>/£' 



dr^ ldx v - \dr^\ \dr u 



(42) 



with = |r M — r u \. For & <C 1 the first term in (42) reproduces the magnetic repulsion of 
London vortices with A' = A/yl — b ~ A; this magnetic interaction is vectorial due to the 
product dr^ ■ dr u containing the cosine of the angle between two line elements. The second term 



of shorter range £' = £/y2(l — b) may be interpreted as an attraction caused by the overlap 
of the vortex cores, the regions where the order parameter is reduced: two overlapping cores 
require less (positive) condensation energy than two separated cores, thus the cores attract. 
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Figure 5: Two profiles of the magnetic field B(x,y) and order parameter \i/j(x, y)\ 2 along the 
x axis (nearest neighbor direction) for triangular vortex lattices with lattice spacings a = 4A 
(b = 0.073, bold lines) and a = 2A (b = 0.018, thin lines). The dashed line shows the magnetic 
field of the isolated flux line from Fig. 4. From Ginzburg-Landau theory for k — 5. 



This attraction has scalar character, hence the product \dr^\ | rfr^ | . The attractive second term 
in (42) removes the logarithmic divergence of the magnetic repulsion at zero distance, since both 
terms have the same singularity but of opposite sign, such that they cancel. The total potential 
(42) is thus a smooth function that at r^ u = starts with a finite value and then decreases 
monotonically to zero with increasing distance — > oo. 

For straight parallel vortex lines the general 3D energy expression (42) simplifies to the sum 
of the vortex self energies, F se if = § B c i/fi per unit length and per vortex, and the interaction 
energy F int of all vortices per unit length, 



2 




(43) 



Here K (x) is a modified Bessel function, see Eq. (21). The effective 2D interaction potential in 
(43) is a smooth, monotonically decreasing function with a finite value at r^ v = since the two 
logarithmic singularities of the K functions cancel each other. As in the 3D expression (42), 
the first term in (43) is the magnetic repulsion of the straight vortices, and the second term is 
an attraction due to gain in condensation energy during the overlap of vortex cores. 



9 Vortex lattice solution for all k> and B 

Abrikosov's solution method for the periodic vortex lattice starts from the linearized GL theory 
and is thus valid only at large inductions B near the upper critical field B c2 . First numerical 
solutions for all B and k were obtained by the "circular cell method" [H] that approximates the 
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Figure 6: Countour lines of u(x,y) = \ip\ 2 and B(x,y) for k = 5, b = 0.5. 



hexagonal Wigner-Seitz cell of the triangular vortex lattice by a circle and solves a cylindrically 
symmetric problem, see also Ref. [15] . The periodic solution in the entire ranges of reduced 
induction < b = B/B C 2 < 1 and GL parameter l/v2 < k < oo may be obtained for 
bulk superconductors by the following numerical method [TSJ, [TBI E] • We start from the free 
energy functional F, Eq. (2), and minimize it with respect to the real and periodic functions 
u>(x,y) = f 2 = \ip\ 2 (order parameter) and Q(x, y) = A — Vip/n (negative super velocity) or 
zB(x,y) = V x Q (induction). We consider periodic lattices with one flux quantum per vortex. 
In the sense of a Ritz variational method we use Fourier series for the periodic trial functions 
with a finite number of Fourier coefficients <2k and 6k, 



w(r) = $>K(l-cosKr), (44 ) 

K 

B(r) = fi + ^&KCOsKr, (45) 



K 



Q(r) = Q^ + ^K^^sinKr, (46) 

where r = (x,y) and K = (K x ,K y ) are the reciprocal lattice vectors (8) of the vortex lattice 
with positions (7). In all sums here and below the term K = is excluded. In (46) Qa{ x , y) is 
the super velocity of the Abrikosov B c2 solution, which satisfies 

VxQ,= [B-$oE J 2(r-R)]z, (47) 

R 

where ^(r) = S(x)8(y) is the 2D delta function. This relation shows that is the velocity 
field of a lattice of ideal vortex lines but with zero average rotation. Close to each vortex center 
one has Qa(i") ~ r' x z/(2/tr /2 ) and ui(r) oc r' 2 with r' = r — R. In principle Qa(i") may be 
expressed as a slowly converging Fourier series by integrating (47) using divQ = divQ^ = as 
in Ref. [TB]. But it is more convenient to take from the exact relation 

2 K 10 a 

where uy\,(x,y) is the Abrikosov B c2 solution given by the rapidly converging series (11). With 
(48) the numerical method becomes highly accurate. Note that the ansatz (46) assumes that 
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Figure 7: The magnetic field variance o = ([B(x, y) — B] 2 ) of the triangular FLL for k = 0.85, 1, 
. . ., 200 plotted in units of B c2 as y/a- (n 2 — 0.069)/£> C 2 (solid lines) such that the curves for all k 
collapse near 6 = 1. The dashed lines show the same functions divided by (1 — 6) such that they 

tend to a finite constant 0.172 at 6 = I. All curves are plotted versus \/b = \J B/B c2 to stretch 
them at small 6 values and show that they go to zero linearly. The upper frame 0.383 is the 
usual London approximation. The limit for very small 6 is shown as two dash-dotted straight 
lines for k = 5 and k = 10. The upper frame 0.383 shows the usual London approximation. 

divQ = 0. This assumption can be shown to be exact at high and low inductions, but I did 
not find a proof that it is true in the general case, though it is satisfied numerically with high 
precision for the periodic vortex lattice at all B and k. 

The solutions ui(r) and B(r) may be computed by using a finite number of Fourier coefficients 
ok and 6k and minimizing the free energy F(B, k, cik, 6k) with respect to these coefficients as 
done in [16J. However, a much faster and more accurate solution method [T5J [17] is to iterate 
the two GL equations 8F/8u = and 8F/5Q, = written in appropriate form. The iteration is 
stable and converges rapidly if one isolates a term (—V 2 + const) (lo, Q) on the l.h.s. and puts 
the remaining terms to the r.h.s. as an "inhomogeneity" of such London-like equations, e.g., 

{-V 2 + 2k 2 )uj = 2k 2 (2u -co 2 - ujQ 2 - g) , (49) 
(-V 2 + cu)Q 6 = -ljQ a - (u -u)Qb, (50) 

with the abbreviations g(r) = (Va>) 2 /(4k 2 cg>), Q b = Q — Q A , V x Q b = B(r) — B, and to = 
(lo) = Y,K a K- Equations (49), (50) introduce some "penetration depths" (2k 2 )" 1 / 2 = ^j\pi 
and u) -1 / 2 = \/u) l l 2 (in real units), which stabilize the convergence of the iteration. Acting on 
the Fourier series u (44) and (46) the Laplacian operator V 2 yields a factor —K 2 , which 
facilitates the inversion of (49) and (50). Using the orthonormality 

2 (cos Kr cos K'r) = 5 K k' (51) 

valid for K / 0, one obtains from (44), (45) cik = — 2(o;(r) cosKr) and 6k = 2(B(r) cos Kr). 
The convergence of the iteration is considerably improved by adding a third equation which 



13 



o 








0.2 



0.4 



0.6 



0.8 



b = B/B 



c2 



Figure 8: The shear modulus of the triangular vortex lattice in bulk superconductors as 
function of the reduced induction b = B/B c2 for GL parameters k = 0.4, 0.5, 0.6, .707, 0.75, 1, 
1.4, 2, 3, 5, 7, 10, 100, in units ^/(lOOO/xo). For « < 2" 1 / 2 = 0.707 one formally has negative 
shear modulus cqq < 0, though vortices and a vortex lattice are energetically not favorable in 
bulk type-I superconductors. 

minimizes F (2) with respect to the amplitude of u, i.e., OF /duo = 0. This step gives the largest 
decrease of F. The resulting three iteration equations for the parameters ax and b K then read 



withp = (Vw x Q)z = Q y duj/dx — Q x duj/dy and g = (Vcj) 2 /(4k 2 cj) = (V/) 2 //t 2 as above. The 
solutions uj(t), B(r), and Q(r) are then obtained by starting, e.g., with ok = (1 — b) a K [the 
Abrikosov solution (11)] and 6k = and then iterating the three equations (52)-(54) by turns 
until the coefficients do not change any more. After typically 25 such triple steps, the solution 
stays constant to all 15 digits and the GL equations are exactly satisfied. 

Since all terms in (52) - (54) are smooth periodic functions of r, high accuracy is achieved by 
using a regular spatial 2D grid, e.g., Xi = {i — \/2)x\/N x (i — 1 . . . N x ) and yj = (j — l/2)y 2 / (2N y ) 
(j = l...N y , 2N y ss N^/xi) with constant weights x 1 /N x and y 2 /(2N y ). These N = N x N y 
= 100 to 5000 grid points fill the rectangular basic area < x < xi, < y < y 2 /2, which is 
valid for any unit cell with the shape of a parallelogram. Spatial averaging (...) then just means 
summing iV terms and dividing by N. Best accuracy is achieved by considering all K mra vectors 
within a half circle |K mn | < K m3iX , with 20N/(x 1 y 2 ) chosen such that the number of the 

K mn is slightly less than the number iV of grid points. The high precision of this method may 
be checked with the identity B(x,y) / B c2 = 1 — uj(x,y), which is valid at k — l/y/2 for all b. 
This relation is confirmed with an error < 10~ 9 . The reversible magnetization M = B — B a and 



4k 2 ((uj 2 + ojQ 2 -2u + g) cos Kr) 
K 2 + 2k 2 



(52) 
(53) 
(54) 



a K • {u - ujQ 2 - g) / (u 2 ) , 
-2{[{u - u)B{r) +p]cosKr) 



K 2 + u 
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Figure 9: Magnetic field lines for a superconductor film calculated from Ginzburg-Landau theory 
for the triangular vortex lattice. Shown is the example b = B/B C 2 = 0.04, k = 1.4, triangular 
lattice with vortex spacing (unit length) Xl = 3- 1 / 4 (2$ / J B) 1 / 2 = 5xi(B c2 ) « 10A, film thickness 
d = 0.8xi ~ 8A. The left half shows the field lines that would apply if the field inside the film 
would not change near the surfaces z = ±d/2 marked by dashed lines. The right half shows the 
correct solution. The density of the depicted field lines is proportional to |B(r)|. 

the equilibrium field B a = n^dF/dB (the applied field) are easily computed from Doria's virial 
theorem [18], which in our reduced units reads 

_ {p-r + 2B{ X) yf) 
Ba ~ 2(B) • (55j 

In this way we find the lower critical field [15J, B c i(k) = limg_^ B a (B, k), 

$ ri , ( S1 , B cl In « + £*(«) 



B c1 (k) = — [ln« + «(«)], h cl ^ 2 



cx(k) = ctoo + exp[— Co — Ci In k — C2(ln k) ] ± e (56) 

with ctoo = 0.49693, c = 0.41477, c x = 0.775, c 2 = 0.1303, and e < 0.00076. This expression 
yields at K — l/v2 the correct value h c \ = 1 and for k ^> 1 it has the limit a = 0.49693. A 
simpler expression for qz(k), yielding an h cl with error still less than 1% and with the correct 
limits at k = l/y/2 and K 3> 1, is 

a(k) =0.5 + (l + ln2)/(2K-V2 + 2). (57) 

The resulting magnetization curves M = B — B a are shown in Fig. 3. They are well fitted by 

h(b, k) = —— w h c \ + 



B c2 l + c 2 6 + c 3 6 2 ' 

Cl = {l-h cl f/{h cl -p)\ 
c 2 = (1 - 3/i c i + 2p)/(h cl -p) , 

c 3 = l + (l-/i cl )(l-2/i cl + p)/(/i cl -p) 2 (5? 
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Figure 10: Profiles of order parameter u>(x,0,z) and magnetic field B z (x,0,z) for the case of 
Fig. 9, film thickness d = 0.8x\ ~ 8A. The solid lines show u and B in the center of the film 
(z = 0) and the dashed lines at the film surfaces. The dotted line indicates the average induction 
B equal to the applied field B a . 

with h cl from Eq. (56) and p = -dm/dfc| 6=1 = l/[(2/t 2 - 1)(3 A + 1] (m = b - h = M/B c2 ), 
(3a = 1.15960 (1.18034) for the triangular (square) vortex lattice. This h(b) satisfies the exact 
relations: h(0) = h cl , h'(0) = h"(0) = h"(l) = 0, h{l) = 1, h'{l) = 1 - p(«). The M(B a ) 
from the fit (58) applies for not too large k < 10 . . .20. For larger k, better fits are given in 
Ref. [15J, where it is also shown that the often used "logarithmic law at B c \ <C B a <C £? c2 " for 
the magnetization M(B a ) = B — B a has very limited range of validity. 

Figure 4 shows the profiles B(r) and u>(r) = \ip{r)\ 2 for an isolated vortex line B — > from 
GL theory for k — 2, 5, and 20. Profiles for the vortex lattice are shown in Fig. 5 for k = 5 
and for two inductions at which the vortex spacing is a = 2 A and a = 4A. Contour lines for 
u(x, y) and B(x, y) are plotted in Fig. 6 for k = 5 at b = 0.5. At b > 0.7 the contours of u(x, y) 
(see right hand part in Fig. 2) and B(x, y) practically coincide, and at b < 0.3 the contours are 
nearly circular, around well separated vortex cores and field peaks. 

The variance of the magnetic field 



is shown in Fig. 7. In the low- field range 0.13/k 2 b 1 one has for the triangular lattice 
the London limit a = 0.00371<I ) q/A 4 (upper frame in Fig. 7), at very small b <C 0.13/k 2 one 
has a = (&fc 2 /87r 2 )$g/A 4 (dash-dotted straight lines in in Fig. 7), and near b = 1 one has the 
Abrikosov limit a = 7.52 ■ 10" 4 ($g/A 4 )[K 2 (l - 6)/(k 2 - 0.069)] 2 . This field variance is needed, 
e.g., for the interpretation of Muon Spin Rotation (/xSR) experiments [T9l 1201 [2T| [22]. 

In Fig. 8 the shear modulus of the bulk triangular vortex lattice is plotted versus the reduced 
induction for various GL parameters k. Note that for k = l/y/2 where B cl = B c = B c2 = 
$ /(47tA 2 ), one has c 66 = 0. One can show that in the particular case k = l/y/2 all possible 
vortex configurations have the same free energy F = BB c \j '//o, e.g., triangular and square lattice, 
lattices with two flux quanta per vortex, or all vortices merged into one giant vortex. 



a = ([B(x,y)-B} 2 )= £ 5 2 



(59) 
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Figure 11: The shear modulus Cqq of the triangular vortex lattice in films with thicknesses 
= 0.1, 0.32, 0.56, 1, 1.8, 3.2, 5.6, 10, and 32, plotted versus b for k = 0.5. This c 66 
is positive, i.e., the triangular vortex lattice is stable, for sufficiently thin films or for small 
inductions. For d ^> £ the bulk cqq at the same k = 0.5 is reached (dash-dotted line), and for 
d <C £ the bulk c 66 in the limit k ^> 1 is reached (dashed line). 

10 Vortex lattice for thin and thick films 



The 2D Fourier method for the bulk vortex lattice in Sec. 9 can be generalized to the 3D problem 
of vortex lattices in infinite films of arbitrary thickness d put into a uniform magnetic field B a , 
since then the functions u>(x, y,z) and B(x, y,z) are still periodic in the (x, y) plane. I consider 
here the case when B a = B a z is perpendicular to the film plane [23], though in principle the 
Fourier method applies also to tilted applied field. The total free energy F tot per unit volume of 
the infinite film is the free energy, Eq. (2), plus the stray-field energy -F stra y, i.e., the energy of 
the magnetic field variations outside the film, 

F tot = F + , F stray = 2 / (B(r) 2 - B 2 ) x , y dz . (60) 

d Jd/2 

One has B = B a since all field lines have to cross the infinite film. The factor of 2 in (60) comes 
from the two half spaces above and below the film, which contribute equally to -F s tra y - The stray 
field B(x, y,z > d/2) with constant planar average (B(x,y, z)) x ^ y = Bz is determined by the 
Laplace equation V 2 B = (since V • B = and V x B = in vacuum) and by its perpendicular 
component at the film surface z = d/2, since B z has to be continuous across this surface. The 
trial functions for u>(r), B(r) = £>z + b(r), and Q(r) = Qa(x, y) + q(r) inside the film (\z\ < d/2) 
are now 3D Fourier series [23J, 

oj{y) = q K (l — cos Kj_r_|_) cos K z z , 

K 

^( r ) = X/ ^ K CQS K-L r -L cos K z z , 
K 
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Figure 12: The magnetization M of infinite films of thickness d/£ = 0.1, 1 ,3, 10, oo with a 
triangular vortex lattice generated by a perpendicular magnetic field B a . Plotted is —M/B C 2 
versus b = B/B c2 = h = B a /B c2 for k = 0.5, 0.707, 1, 1.5. 

K K ± 

C l( r ) = H b K Z ^2~^ Sil1 K ^ r -L cos K zZ , (61) 
K K ± 

Here r = (x,y,z), r± = (x,y), K = (K x ,K y ,K z ) with K ± = (K x ,K y ) from Eq. (8), and lf 2 = 
(2ir/d)l, I = 0, 1, 2, . . .. In all sums here and below the term Ki = is excluded. Minimizing F tot 
with respect to the coefficients ok and 6k and using the appropriate orthonormality relations 
one arrives at iteration equations for the ok and 6k similar to Eqs. (52)-(54). The solution is 
then obtained by first finding the 2D bulk solution as in Sec. 10 by considering only the terms 
with K z — 0. The magnetic field lines then still have an unphysical sharp bend at the surface, 
see left part of Fig. 9. Next we allow for the terms with K z ^ 0. This yields a "mushrooming" of 
the field lines of each vortex when it approaches the film surface such that these lines smoothly 
cross the surface with no bend, see right part of Fig. 9. 

Profiles u(x,0,z) and B z (x, 0, z) are shown in Fig. 10 for z = (middle plane of the film) 
and z — d/2 (film surface). One can see that the spatial variation of B at the surface is reduced 
from its bulk value by nearly 1/2. Outside the film, the transverse field components B x , B y 
rapidly decrease as exp(— |K 10 |z') w exp(— 2nz' /a) where z' = \z\ — d/2 is the distance from the 
surface. Interestingly, the profile of the order parameter uo in films is almost independent of z. 

The shear modulus cqq of the triangular vortex lattice in films can be positive even when 
k < l/y/2, provided the film thickness d is smaller than the coherence length £, see Fig. 11. This 
means that a stable vortex lattice may exist in thin type-I superconductor films. Our numerical 
result confirms the Cqq of films that was calculated analytically near B c2 by Albert Schmid 



In such infinitely extended films one has B = B a since all field lines have to pass the film. 
Therefore, the magnetization M of the film cannot be calculated as a difference of fields, but one 
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has to take the derivative of the total free energy, M = B — dF tot /dB. A more elegant method 
calculates M by Doria's virial theorem [18], which for bulk superconductors yields Eq. (55). 
Indeed, it was shown recently [25] that this virial idea can be generalized to films of arbitrary 
thickness and M may be calculated directly from the GL solution for the film, with no need to 
take an energy derivative. The resulting magnetization of the film is plotted versus the applied 
field in Fig. 12 for various film thicknesses d and various GL parameters k. The curves for various 
d cut each other at b « 0.1/ k. For k = 0.707 the thick film limit (d ^> £ but still film width 
w ^> d) yields a straight line, —M/B C 2 = 1 — 6. This result is valid for large demagnetization 
factor N, 1 — N <C 1; it differs from the bulk result for M shown in Fig. 3, which is valid for the 
demagnetization-free limit of N <C 1. 

11 Remarks on the magnetization 

One may ask why the magnetization M = m/V is not calculated via the general definition 
of the magnetic moment m = V~NL = | J v d 3 rr x j [221 EZ] from the current density j(r), 
which is easily calculated as a periodic function by our Fourier method, both for bulk and film 
superconductors. However, into this definition enters not only the periodic part, i.e., the vortex 
currents circulating inside each vortex cell; this contribution even would give the wrong sign of 
m. The main contribution to the magnetic moment of a superconductor of any shape comes 
from the screening currents that flow near the surface of the specimen. The magnetization in 
superconductors is thus not a volume property as it is in magnets. For example, for a long 
cylinder in parallel field B a the magnetization M = B — B a < is composed of the positive 
contribution B of the vortex currents and the negative (and larger) contribution B a of the 
surface currents that screen the cylinder from the applied field B a before vortices are allowed 
to penetrate. Near B c2 , both terms nearly compensate and \M\ is a small difference of two big 
terms. In the film geometry, and for most other shapes of the superconductor, the screening 
current is not easily known but has to be computed; such computations for thin and thick strips, 
disks, and plates of macroscopic size 3> £ are presented in [26J. Landau knew this problem, since 
he worked on demagnetization factors and on the intermediate state in type-I superconductors 
[27] , and he used the thermodynamic definition of the magnetization as an energy derivative. In 
my view, the simple formula (55) derived by Doria, Gubernatis and Rainer [IB] by scaling the 
coordinates in GL theory and finding a novel virial relationship between kinetic and potential 
energies from which the equilibrium field B a = B — M follows, was a fundamental discovery that 
occurred long after the publication of GL theory [TJ . 
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